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■rj" ' Abstract 

' The evaluation of the path-integral representation for stochastic processes in the 

. weak-noise limit shows that these systems are governed by a set of equations which are 

2 ' those of a classical dynamics. We show that, even when the noise is colored, these may 

- ^ , be put into a Hamiltonian form which leads to better insights and improved numeri- 

■ cal treatments. We concentrate on solving Hamilton's equations over an infinite time 

interval, in order to determine the leading order contribution to the mean escape time 
for a bistable potential. The paths may be oscillatory and inherently unstable, in which 
case one must use a multiple shooting numerical technique over a truncated time period 
in order to calculate the infinite time optimal paths to a given accuracy. We look at 
two systems in some detail: the underdamped Langevin equation driven by external ex- 
ponentially correlated noise and the overdamped Langevin equation driven by external 
quasi-monochromatic noise. We deduce that caustics, focusing and bifurcation of the 
optimal path are general features of all but the simplest stochastic processes. 



1 Introduction 

The subject of noise induced activation has received a great deal of attention in the last 
decade or so with the development of new techniques which allow systems where the noise 
is not white (i.e. colored) to be studied in a systematic and controlled way While 
these techniques were being refined it was natural that only the simplest of systems were 
studied: those which consisted of a single particle moving in a one-dimensional potential 
with a relatively simple form of noise. More recently the investigation of models acted upon 
by white noise and having more than one degree of freedom has revealed novel effects such 
as caustics and focusing singularities 0-0. In these systems it was also found that the 
leading order term in the expression for the mean escape time, that is the action, was reduced 
unexpectedly 0. 

The appearance of such features can be understood in the following way. In the limit of 
weak noise the dynamics of the system is governed by a set of equations which are the extrema 
of the action in the path-integral formulation of the stochastic process. These equations have 
the same form as those of Newtonian mechanics, if the original stochastic dynamics was 
underdamped and the noise was white. Hence the dynamics of these stochastic processes is 
controlled by trajectories in a 2n-dimensional phase space, where n is the number of degrees 
of freedom of the system. Placed in this context, it is not unnatural to expect caustics and 
focusing singularities [0. In this paper we will show that similar effects will also occur in 
systems with one degree of freedom, but with a more complicated type of noise. We do this 
by using a generalized Hamiltonian formalism to show that a phase space can be constructed 
which is multi- dimensional and hence, since these systems do not satisfy detailed balance, 
will be expected to show the phenomena mentioned above. 
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Our starting point is the observation that a system consisting of a single degree of freedom, 
but acted upon by a rather general form of external noise 0, can be written as a Markov 
process which consists of a number of equations, only one of which involves a noise term 
(which is white). These equations can be combined into a single equation, at the expense of 
introducing higher time derivatives. Essentially we have traded a simple system acted upon 
by a complicated noise term, for a complicated system acted upon by a simple white noise. We 
shall give explicit examples later in this paper. A process of the kind we have been describing 
can be defined by the generic stochastic differential equation 

f{x,x,x,...,x^-^;t)=r]it) (1) 

where 77 (t) is Gaussian white noise of strength D. 

Let us now make the above comments on the emergence of a classical dynamics in the weak 
noise limit more concrete by outlining how a path-integral representation for the conditional 
probability distribution of a process defined by (0) can be written down ||^. One begins by 
using (|l]) to transform the probability density functional for white noise given by 

P[r^]=Cexp[-^j\\t)dt^ (2) 

to the probability density functional for that of the coordinate x: 

P[x] = ^J[x] exp (^-^ j (3) 

where S[x] is the action mentioned above and J[x] is the Jacobian of the transformation. The 
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action is so called since it may be written as 



S[x] = dt L{x, x,x, . . . , 



X 




(4) 



where L(x, . . .) given by 



, . . . , X^ \ t) ^ [fi^i X, X, . . . ,x^ \ 



(5) 



This has the form of a Lagrangian for a mechanical system, if the noise is white and the 
motion is overdamped, so that no time-derivatives higher than the first appear in (|l|). For 
colored noise processes of the type we are investigating here, there are higher time-derivatives 
in the Lagrangian and the analogy is now with a generalized form of mechanics. The precise 
form of the Jacobian factor will not be required in this paper since we will be performing our 
calculations to leading order only, and the Jacobian only enters at next order. Probability 
distributions, correlation functions and other quantities of interest can be found by integration 
of the appropriate functions over paths x{t) with weight (P). In the limit of D ^ these path- 
integrals can be evaluated by the method of steepest descents; the paths which dominate the 
integrals being the ones for which 6S[x]/6x = 0. This leads to the Euler-Lagrange equation 
for the optimum path which will be in general a 2n^^ order non-linear differential equation 
given by 



In general, this equation will have no analytical solution and one has to rely on numerical 



2n coupled first order non-linear differential equations. It would be convenient to derive the an 




(6) 



techniques |TU[ . A numerical solution will involve the decomposition of such an equation into 
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expression for the optimal path in such a format automatically. This is instantly provided by 
using Hamilton's formalism as an alternative to the Lagrangian method. Another advantage 
of using this method is that the greater geometrical structure of Hamiltonian mechanics gives 
one a better insight into why the optimum paths take on the particular form that they do. 

The outline of the paper is as follows. In section two we construct the generalized Hamil- 
tonian formalism appropriate to problems of this type and give the case of the underdamped 
Langevin equation driven by white noise as an example. In section 3 we use the formalism to 
find the mean first passage time for this underdamped problem, but now with exponentially 



correlated noise; a task which could not be achieved using the Lagrangian formalism |10 



The case of quasi-monochromatic noise is discussed in section 4 and we conclude in section 5. 

2 Hamiltonian Formalism 

For a dynamical system which is defined by a Lagrangian of the form L(x, x,x, . . . , x*^"'-'; t), a 
Hamiltonian structure can still be constructed (see, for instance, |TT]]). To do so one intro- 
duces a generalized coordinate vector q spanning an n dimensional space with components 
{gi, . . . , Qn} such that 

q^ = x^'-'^ (7) 

and one writes 

n 

L(a;,i;,x, ...,x(");t) =J2piqi-H{q,p;t) (8) 

i=l 

where the pi's have yet to be defined. Now if one demands that pi = —dH/dqi, it follows 
from (D that 

dL 

Pj+Pj+i ;j = 0, ...,n (9) 



where po and Pn+i are defined to be zero. From (H) and (j^) one sees that the components 
{pi, . . . ,pn} of the generalized momentum vector p should be taken to be 

p. - p-^y-'§^. {^) (1°) 

Hence, by construction, the optimum path given by the 2n*^ order differential equation @ 
can also be found by solving the 2n first order differential equations 



(U) 

dpi 



and 

dH , , 

= 12 

dqi 

If the Lagrangian does not involve time explicitly, then the Hamiltonian, H, also has no 
explicit time dependence, and since dH/dt = dH/dt, the Hamiltonian is conserved. This 
reduces by one the number of integrals that have to be performed. 

As an example we shall consider the underdamped Langevin equation driven by white 
noise 

mx + ax + V'{x) = rjit) (13) 

where r]{t) is Gaussian white noise of strength D. Here V{x) is assumed to be a double well 
potential and a is a friction constant which will be set equal to unity by an appropriate choice 
of units of time. The Lagrangian for this process is given by 



L{x,x,x) = ^{mx + X + V'{x)f (14) 
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and the Hamiltonian is found to be 



Hici, p) = p,q, + 4 - - (g2 + V'iq,)) (15) 



The optimum path is then the solution of Hamilton's equations: 



qi = q2 



?2 



2P2 q2 V'iqi] 



m? m m 



. P2V"{qi) 

Pi = 

m 



P2 = Pi (16) 



m 



and the action is given by 



S= f ^dt (17) 



to m? 



In this case we can find the required solutions explicitly enough to allow us to write down the 
action in closed form. We are searching for solutions which begin at extrema of the potential 
with all time derivatives of the coordinate equal to zero. This immediately tells us that H = 
for these solutions, which is a common feature in models of this type. For this simple case 
there are only two of these solutions: a "downhill" solution given by mx + x + V'{x) =0 and 



an "uphill" solution given by mx — x + V'{x) = [10]. These solutions can easily be found as 



H = solutions of (p^Gf): the downhill solution has Pi = P2 = and zero action and the uphill 
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solution has pi = V'{qi), p2 = mq2 and action given by 



1 1 



(18) 



The interpretation of (|1^) depends on exactly what quantity is being calculated. For example 
if one wished to find the stationary probability distribution then one would take to ^ —oo 
so that 5* = mx'^/2 + V{x) in terms of the original variable x, which just gives the Maxwell- 
Boltzmann distribution. On the other hand, if one wished to find the mean escape rate from 
a potential well, one is interested in paths which take an infinite time to interpolate between 
stable and unstable points of the potential and are at rest at both ends. This gives S = AV, 
where AV^ is the barrier height. For the rest of this paper we will restrict ourselves to the 
calculation of this quantity and so will assume an infinite time interval in what follows. Since 
for colored noise processes, which are the real interest of this paper, we cannot, in general, 
calculate the action explicitly, we will choose the specific double-well potential 



to illustrate our techniques. If we choose to investigate activation from the left-hand well to 
the right-hand well, then the section from zero to plus one will be a downhill path with zero 
action. Thus we need only concern ourselves with the section of the path from minus one 
to zero. Having illustrated the technique on a simple white noise problem we now go on to 
investigate the same system, but acted upon by exponentially correlated noise. 




(19) 
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3 Exponentially correlated noise 

In this section we consider the process modelled by the Langevin equation 



mi + x + V'{x) = ^{t) (20) 



where ^(t) is Gaussian colored noise whose correlation function is given by 



mm) = ^^wP-^) (21) 



This represents the simplest generalization of the noise in the system first investigated by 



Kramers [12], which it reduces to in the r ^ limit. One can see that it is the simplest 
generalization by replacing ( pTj) by the condition that ^ obeys the first order differential 
equation 

ri + ^ = v{t) (22) 

where is a Gaussian white noise of strength D. Equations (^) and (^) form an equivalent 
Markov process with two degrees of freedom. These equations may be combined into the 
single third-order stochastic differential equation 

mx + X + V'{x) + T{mx + x + xV"{x)) = rj{t) (23) 



This is of the form (0) with n = 3 and so we expect to be able to describe the weak-noise 
limit of this system using either Lagrangian or Hamiltonian dynamics. 

The Lagrangian approach to this problem has been investigated by Newman et al |TD 



However these authors were only able to explore the dynamics of the system for relatively 



small masses; they were unable to analyse the underdamped regime. We shall show in this 
section that the Hamiltonian approach allows us to do this. 

Using and p^ ) we can write down a Lagrangian for this system: 



1 2 
L(x, X, X, x) = - [(mx + X + V'{x)) + T{mx + x + xV"{x))] (24) 



The equivalent Hamiltonian is found in the way described in section 2 to be: 



2 

H{q^P) = Piq2 + P2qz + ^ - — (mgs + ^2 + V\q,)) - ^{q^ + q2V"{qi)) (25) 



with Hamilton's equations given by (11) and ([T^). The action reduces to 



2 

S = I dt^ (26) 



We now wish to find the value of the infinite time action for the bistable potential (|I^). The 
downhill solution is, as usual, trivial: p„ = 0; n = 1, 2, 3, which gives zero action. The uphill 
solution cannot be found analytically for general m and r and only perturbative methods and 
numerical solutions are available. The m = 0, general r problem is extensively discussed in 
, along with perturbative expansions in the small r and large r regimes. Therefore we will 
restrict ourselves to m > 0. In Ref a numerical calculation of the action for certain values 



of m and r have been given, as well as perturbation expansions for small m and small r. In 
the rest of the section we will expand on this treatment, extending it and investigating the 
previously unexplored underdamped regime. 
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For general m, but small r, the action for the uphill path has the form 



S{m, t) = So + T^S^im) + O(r^) (27) 



where 5*0 = 1/4 is the white noise action for this potential. The first correction Si{m) has 



the simple form |T0 



oo 



Si{m) = ildt (28) 



oo 



where Xq is the optimal path for white noise (r = 0) and is given by the non-linear differential 
equation 

mxQ - xo + V\xq) = (29) 
with the boundary conditions Xo{—oo) = —1 and Xo{oo) = 0. 



In the paper by Newman et al |T0[ this quantity was calculated for small m only. However 
equation ( p9D is stable if integrated backwards in time, i.e. starting at xq = going to 
xq = —1 (stability is discussed further when the full solution is considered). Hence, one can 
use a simple initial value integrating scheme, such as a fourth-order Runga Kutta, starting 
with an infinitely small velocity xq = S. While in a formal sense the path over the infinite time 
interval is only found in the limit 5 0, in practice we find that if S is small, the value of the 
action does not depend on it. The results for Si{m)/So are plotted in figure 1 as a function 
of logiQ^m). The dotted line shows a seventh-order perturbative calculation of 5*1 (m): 



^ , , ^ / 1 m rn? 2rn? 3m'^ 9m^ 3m^ 778m'^\ ^, o, 
S.(m)=S„ - + ---- — + — + — - — -— + 0(m») (30) 



Figure 1 shows the excellent agreement between the series (pOf) and the numerical solution 
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for m less than about 0.3, and the catastrophic failure of the series above that value. This 
breakdown of perturbation theory may be due to the change in the nature of the solutions 
that occurs at m = 1/8 (see below) and a calculation of more terms in the series ( pO| ) might 
show that the value of m at which the breakdown occurs approaches the value 0.125. This 
figure also shows that the value of Si{m) has a maximum when plotted against log^o('^)- 
Such maxima are also seen when plotting actions against log^o(^) (EH Figure 2 below). 
The existence of these maxima are a consequence of the non-linear nature of the problem. 

Now let us go on to a numerical study of the solution of Hamilton's equations for general 
m and r. As a first step we linearize the equations about the endpoints qi = a, where 
a is zero or minus one. To do this we approximate the potential by a parabola V{qi) = 
V{a) + \V"{a){qi — a)^, which leads to linear Hamilton's equations with solutions of the form 
qi = a + J2n ^nC^"* where An are arbitrary constants specified by the boundary conditions. 
The An have six possible values: 

1 / 1 Jl-4mV"{a)\ 

Xn = ±- , ± — ± — (31) 

r \2m 2m I 

When a = 0, V"{a) = —1 and we require that gi — as t — oo for an uphill path, so we 
select only those A„ which are negative. Conversely, when a = — 1, V'(a) = 2 and we require 
that gi — 1 as t ^ — oo, thus we must only take values of A„ which are positive, i.e. 

= a + ^IHl™) , (32) 

' ' r \2m 2m J \2m 2m J 

If m < 1/8 no problem arises. However if m > 1/8, two of the quantities in ( |32D are complex 
which is a signal that the solutions may be oscillatory. Actually if r > 2m the real one 
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dominates, implying that the solutions are not oscillatory near to the stable fixed point in 
this case. In summary, we can say that, as in the underdamped Langevin equation with 
white noise, the system oscillates about the bottom of the potential wells before making a 
transition, unless m < l/8orm> 1/8 and r > 2m. In this case the substitution y{x) = x{t), 
which was the basis of the approach in [^, fails and another technique has to be used. 

The advantage of solving for y{x) is that since — 1 < x < 0, the differential equation has 
to be solved in a finite range. Unfortunately, in the region of parameter space where the 
oscillatory solution exists we have to solve Hamilton's equations over the range —oo<t<oo. 
In practice, of course, we have to truncate this span to a large, but finite, value T and use 
the boundary conditions 

?i(-|) = -l ?i(|) = 

?2(±|)=0 

g3(±|)=0 (33) 

and calculate the action 

The extrapolation T ^ cxd is, in fact, not a problem; the actual transition happens over a 
time scale of order a few m and for the majority of the time the particle is almost at rest 
at the two endpoints. Furthermore, this decay of the position, velocity, etc, at the endpoints 
is exponential, which means that the truncation can be carried out extremely accurately. 
However the simple initial value techniques of solution which are used in shooting routines 
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can no longer be used, as this problem is inherently unstable. This is because, as one can 
see from (|3T|), there are growing solutions at both the end points. In analytic treatments 
these solutions can be ignored by setting the arbitrary constants (An) to zero. Numerically, 
roundup errors introduced either through machine precision or through the solution algorithm, 
make these constants small, but non-zero. Since these solutions grow exponentially, whereas 
the required solution decays exponentially, they soon take over, and any hope of solving the 
problem numerically by this method is destroyed. 

Normally, it is possible to solve inherently unstable problems by using a relaxation tech- 



nique or collocation such as COLSYS |jT3[- However, as noticed in Ref [|T0| this technique 
is poorly convergent when the solution is oscillatory. Instead, one can attempt to proceed 
using either invariant embedding or multiple shooting techniques. It is the latter that we 



have used; calculating the action for r < 2m using MUSN WM. These techniques damp out 



the exponentially growing solution by splitting the total time span into several smaller time 
segments and then matching the solution continuously In this regime it turns out that 
the time of integration, T, needs only to be of the order of a few m in order to obtain reliable 
results. On the other hand, if we try to find the solution as a function x{t) in the regime 
T > 2m, the time of integration needs to be of the order of a few r, which, since we are 
interested in large values of r, becomes a problem. Fortunately, as we have seen, the solution 
can be found as a function y{x) in this case. 

The results of the numerical solution are shown in Figure 2 for several values of the mass. 
It is convenient not to plot the action S{m, r) itself, but the reduced quantity 



S.Kr) = %li^ (35) 
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since this is finite in the hmits r — and r oo, having the values zero and one respectively. 
Here 5*0 is the action when r = and Soo{t) is the action in the large r limit and is given by 

i 

Soo = ^ (36) 

5*0 and 6*00 are m independent. For r > 2m, the hollow points have been calculated using 
{y,x) variables and the solid points using {x,t) variables. The lines are curves of best fit 
through these points to aid the eye. 

This figure shows that numerically the mass has little effect on the action, and the over- 
damped Langevin equation provides an excellent approximation to the action for the under- 
damped system. There seems no reason to expect this a priori, except for the fact that since 
in the limit of small and large r the action is m independent, there is very little freedom at 
intermediate values of r to have significant deviations from the m = result. 



4 Quasi- mono chromatic noise 

In the last two sections two explicit types of noise have been considered: white noise whose 
power spectrum is flat and exponentially correlated noise which has a spectrum centered about 
zero. A type of noise that has a definite color, in the sense that it has a power spectrum peaked 
at a non-zero frequency, is quasi-monochromatic noise (QMN) [|l^-pO|. Systems acted upon 
by QMN are the subject of this section. 
The noise ^(t) is defined by 

^■ + 2re + ^o'^ = r/ (37) 
where rj is a Gaussian white noise of strength D. Hence for the overdamped system x+V'{x) = 

14 



^(t), the Lagrangian is given by 



2V 1 

(x + V\x)) + — (x + xV"{x)) + — (x + xV"(x) 



n 2 



xV'"(x)) 



(3^ 



and the Hamiltonian found to be: 



H{q, p) = p,q, + p2q3 + u^^pI ' Ps {^0(^2 + V) + 2r{q, + q,V") + q,V" + qlV'"] (39) 



The dynamics will be governed, in the limit of weak noise, by solutions of Hamilton's equations 
given by (pJ]) and (|T^). For concreteness we will again consider the potential to be of the 



form (19), and hence we will require the truncated infinite time boundary conditions (^) for 
the uphill path. The action for the uphill path will be given by 



S 



UqpI dt 



(40) 



The downhill path again leads to zero action and will not be considered further. 

An analysis of the linearised Hamilton's equations near the end-points along the lines 
described in the last section, again shows there to be oscillations depending on the value of 
r (in fact oscillations occur for F < min(2,Ci;o)) and so once again we are unable to use the 
(?/, x) parameterization of the solution. Figure 3 shows the generalized coordinates found 
by solving Hamilton's equations for F = 0.45 and cjq = 10. This particular value of ujq 



was chosen to allow comparison with earlier work ||2^ where an approximate solution to the 
classical dynamics was used to calculate the action. The value of F is chosen for clarity: for 
smaller values it is harder to illustrate graphically a complete transition showing the smaller 
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scale oscillatory features characteristic of QMN, whereas for larger values these oscillations 
are absent. The paths have three distinctive features: 

1. An underlying oscillatory factor of angular frequency c^o- 

2. An underlying growth and decay either side of the transition time to given approximately 
by exp(-r|t - to\). 

3. They pass over the top of the potential barrier many times before coming to rest. 

These features only occur if F is less than a critical value Tc (which has a value just less 
than a half), otherwise the solution is that of the system acted upon by white noise to an 
accuracy of 1% (i.e. of order l/uj^). They explain why the approximate treatment given in 



| 20| was successful: there it was assumed that the paths had exactly the features 1) and 2) 
above. The last point mentioned above shows that one has to distinguish clearly between a 
mean first passage and a well transition. 

Figures 4 and 5 show the second and third figures of redrawn with the action calculated 
from the Hamiltonian technique shown as a dotted line. The asterisk on the action S indicates 
that it is the most probable escape path (MPEP) — for F < Fc the escape path can be either 
white-noise-like or oscillatory, but it is the latter that occurs in practice since it has the least 
action and so is most probable. These two figures show the remarkably good agreement 



between solving the full equations and the approximation used in [20|: the value of F^ is 



approximately the same and a maximum value of 5'*/F occurs at F ~ 0.1. From figure 4 one 
can also see that for F < Fc, 5** ~ |F and for F > Fc, S** ~ |. The intersection of these lines 
gives Fc = |, which is a reasonable estimate. 

Difficulties arise for small F as the time required for transition goes as F~^ and hence 
a longer time span, T, is required. If we attempt to rescale time by F, the frequency of 

16 



oscillations now goes as ujq/T, which means we need a finer grid of shooting points to calculate 
the action to sufficient accuracy. So far we have only been able to extend our method down 
to r = 0.05. 

If one writes down an equation for the optimum path (given by equation (12) of pO| ) 
perturbatively in powers of I/cUq' o'^^ ^^"^ ihaX the uphill solution is 

x = V\x) + ol^^ (41) 

which has the corresponding action 

S = \-^ + (—\ (42) 

This approximate solution is independent of the value of F, and exists independently of the 
value of r. For F > Fc it is the global minimum, however for F < F^ one finds that the 
optimum path from minus one to zero bifurcates. In this case this white-noise-type path 
becomes a local maximum and the oscillatory- type path becomes a local minimum. The 
existence of these latter paths is not obvious when solving Hamilton's equations numerically; 
a very thorough search in phase space is required to find them. This situation is common in 
a system such as this with several degrees of freedom: the existence of caustics and focusing 
gives rise to bifurcations in optimal paths 0, which makes the prediction of the correct action 
difficult. 

An obvious extension of this work is investigate the driving of the underdamped Langevin 
equation ( pOD by QMN ^(t) given by (0). We might expect that driving an equation such as 
this with harmonic noise such as QMN, we would find a problem which is inherently unstable 
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with oscillatory solutions. This is indeed the case, and the problem has to be solved using 
the time as independent variable and by use of a multiple shooting technique. Another added 
problem is the introduction of two more coordinates, since the Lagrangian now has fourth- 
order time derivatives. Though this does not cause any further instabilities, it does add to 
the complexity of the problem and further complicates finding the required solution. We shall 
not pursue this extension any further, since it does not introduce any novel features. 

5 Conclusion 

The Hamiltonian formahsm has proved effective for obtaining results for stochastic systems 
governed by complicated differential equations. It has allowed us to understand why optimum 
paths take on particular forms. It also indicates that caustics and focusing appear as general 
features of systems governed by colored noise (even those with only one degree of freedom) 
and not just those with white noise and more than one degree of freedom. The Hamiltonian 
formalism is the natural one in which to investigate and understand these caustics system- 
atically. We have also shown that the technique of of multiple shooting, though slower than 
relaxation and less convergent, has allowed us to study regions which have so far been elusive 
and has opened up the solution of these instanton paths in terms of the original {x, t) vari- 
ables. This could be useful when investigating time-dependent Lagrangians or more complex 
oscillatory problems. We now feel that the structure and general features of optimal paths are 
better understood, and that as a consequence the weak-noise evaluation of escape rates and 
the stationary probability distribution for many stochastic processes is now becoming more 
straightforward. 
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Figure Captions 

1. Si{m) plotted against logiQ(m) 

2. Sr{m,T) against logio(T) for different values of m 

3. qi,q2 and against i for F = 0.45 and ujq = 10 

4. Minimum QMN action S* against F. The dotted line is from the Hamiltonian method 

5. Minimum QMN action S* /T against F. The dotted line is from the Hamiltonian method 
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